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Abstract 



We study the phase diagram and critical behavior of the one- 
dimensional pair contact process (PCP) with a particle source using 
cluster approximations and extensive simulations. The source creates 
isolated particles only, not pairs, and so couples not to the order pa- 
rameter (the pair density) but to a non-ordering field, whose state in- 
fluences the evolution of the order parameter. While the critical point 
Pc shows a singular dependence on the source intensity, the critical ex- 
ponents appear to be unaffected by the presence of the source, except 
possibly for a small change in /3. In the course of our study we obtain 
high-precision values for the critical exponents of the standard PCP, 
confirming directed-percolation-like scaling. 
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I. INTRODUCTION 



Critical phenomena at absorbing-state phase transitions are of longstanding in- 
terest in statistical physics [1-3], being found, for example, in models of epidemics 
[4], catalytic kinetics [5,6], surface growth [7], self-organized criticality [8-11], and 
turbulence [12-14]. The transition between active and absorbing states arises out 
of a conflict between two opposing processes (e.g., creation and annihilation); when 
continuous (as is often the case), it falls generically in the DP universality class 
[15-17]. When two or more absorbing states exist and are connected by a symmetry 
operation, as in branching and annihilating random walks, a new kind of critical 
behavior appears [18-21]. 

Unusual critical behavior also arises in models that can become trapped in one of 
an infinite number of absorbing configurations (IN AC). (More precisely, the number 
of absorbing configurations grows exponentially with the system size. There is no 
special symmetry linking the different absorbing configurations.) Models of this 
sort were introduced to describe surface catalysis [22,23]; their critical properties 
have been studied in detail by various workers [24-31]. In one dimension, the pair 
contact process (PCP) [25], and other models with INAC exhibit static critical 
behavior in the DP class [26,32], but the critical exponents S and 77, associated with 
the spread of activity from a localized seed, vary continuously with the particle 
density in the environment [26,27,31,33-35]. [These exponents are defined via 
the asymptotic (t 00) power laws: survival probability P(t) ~ t~^, and mean 
activity n{t) ~ f^; note also that the mean-square distance of activity from the seed, 
R'^{t) ~ t^.] This anomalous aspect of critical spreading for INAC can be traced to 
a long memory in the dynamics of the order parameter, p, arising from a coupling 
to an auxiliary field (the local particle density, 0), that remains frozen in regions 
where p = [28,30,31]. Theoretical understanding of models with INAC remains 
incomplete. Mufioz ct al. proposed a continuum theory for models with INAC, and 
showed that the static critical behavior is that of DP [28] . Formally eliminating the 
auxiliary field, they obtained a closed equation for the order parameter, in which a 
memory term appears; simulations of this theory also show nonuniversal exponents 
[33]. While a theoretical analysis of variable spreading exponents in the PCP is 
lacking, an analogous phenomenon has been analyzed exactly in the simpler cases 
of a random walk with a movable partial reflector [36] , or of compact DP confined 
to a parabola [37] . 

In previous studies of the PCP, and, indeed, of all nonequilibrium models pos- 
sessing an auxiliary or non-ordering fleld 0, the latter has been allowed to relax to its 
stationary value [38]. Only the effect of varying in the initial condition has been 
explored. Given the surprising results of these studies, it is of interest to investigate 
the consequences of changing the non-ordering field in the stationary state as well. 
To this end, we introduce an external field h that couples to 0, but not to the order 
parameter itself. In concrete terms, h is the rate (per site) of attempted insertions 
of isolated particles. The source may only insert a particle at a vacant site, both of 
whose neighbors are also vacant; in this way no pairs are created, and the absorbing 
nature of configurations devoid of pairs is maintained. In this work we examine the 
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effect of this perturbation of the phase diagram and critical behavior of the PCP, 
using cluster approximations and extensive simulations. 

The balance of this paper is devoted to defining the model (Sec. II); a discussion 
of cluster approximations (Sec. Ill); and analysis of simulation results (Sec. IV). 
We close in Sec. V with a discussion of our findings. 

II. MODEL 

In Jensen's pair contact process (PCP) [25], each site of the one-dimensional 
lattice Z is either vacant or occupied by a particle. Each nearest-neighbor (NN) 
pair of particles has a rate p of mutual annihilation, and a rate 1— p of attempted 
creation. In a creation event involving particles at sites i and a new particle may 
appear (with equal likelihood) at site i—V or at i-|-2, provided the chosen site is vacant. 
(Attempts to place a new particle at an occupied site fail.) In an annihilation event 
a NN pair of particles is removed. The PCP exhibits an active phase for p < pc; 
for p > Pc the system eventually falls into an absorbing configuration devoid of NN 
pairs, but that typically contains a substantial density, 0, of isolated particles. The 
best estimate for the critical parameter is = 0.077090(5) [32]. (Here and in what 
follows, numbers in parentheses denote uncertainties in the last figure or figures.) 

To the above dynamics we now add a third process, addition of isolated particles. 
Each site is bombarded by particles at rate h. (This is a fiuctuating source; the mean 
time between successive addition attempts at a given site is l/h.) An attempt to 
place a particle at site i is successful if and only if sites i — i, and i + 1 are all 
vacant. In the absence of the source, an empty lattice is absorbing. For any h > 0, 
however, the insertion of isolated particles onto an initially empty lattice corresponds 
to random sequential adsorption (RSA) of dimers in one dimension. The saturation 
density for this process is {l-e-'^)/2 = 0.432332... [39]. 

Previous studies leave httle doubt that the static critical behavior of the PCP 
(without a source) belongs to the universality class of directed percolation. Jensen 
and Dickman found that the critical exponents f3 and 7 (which govern, respec- 
tively, the stationary mean of the order parameter, and its variance), and i^w and i'± 
(which govern the divergence of the correlation time and correlation length as one 
approaches the critical point), are all consistent with DP values [26]. In addition to 
these static properties, the exponent 6, which governs the initial decay of the order 
parameter (p oc t^^), starting from a fully occupied lattice, was found consistent 
with DP. More recently, the order-parameter moment ratios were found to be the 
same as those of other models belonging to the DP universality class [32]. 

Starting from a spatially homogeneous distribution of NN pairs (for example, a 
completely filled lattice), the system relaxes to a stationary state. If, by contrast, 
the activity is initially localized (e.g., a single "seed" pair in an otherwise absorbing 
configuration), we may study the spread of activity. As noted above, the critical 
exponents S and 77 characterizing spreading vary continuously with 0, and assume 
DP values only for (^ — (pnat — 0.242 [40,41]. The natural density (f)nat is the mean 
particle density in absorbing configurations generated by the process itself, at Pc, 
starting with all sites occupied. One may, equivalently, define (pnat as the particle 
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density in the critical stationary state, in the thermodynamic hmit. An environment 
with (f) > (pnat favors spreading (and vice- versa), since the higher the particle density 
in the environment, the more pairs will be formed per creation event. 

A kind of spreading phenomenon also arises in the stationary state due to spon- 
taneous fluctuations. In the critical stationary state, we can expect to find inactive 
regions of all sizes; the particle density in large inactive regions is ^nat- When activ- 
ity spreads into such regions, it should follow the same scaling behavior as critical 
spreading with (pnat- Since the exponents governing survival and growth of activity 
in the stationary state are subject to the scaling relations 5 — /3/z/|| and z — 2v^/u\\, 
with /3, and v\\ taking DP values in the stationary state, it follows that the 
spreading exponents take their usual DP values as well, for = 0nat- (In an envi- 
ronment with 7^ 0nat, the advance of activity is no longer equivalent to that in 
the stationary state, and the spreading exponents are not constrained to take DP 
values.) In this work we study static critical behavior with cj) ^ (p^at due to the 
action of the source. 

III. CLUSTER APPROXIMATIONS 

In this section we develop dynamic cluster approximations for the PCP with 
a source of particles. Such approximations often yield qualitatively correct phase 
diagrams [2]. For the standard PCP (no source) the two-site approximation was 
presented by Carlon et al. [42], while the three-site approximation was derived by 
Marques at al. [41]; our results for h = are consistent with these studies. The n-site 
approximation consists of a set of coupled differential equations for the probabilities 
P^^^ of each n-sitc configuration, C. (There are 2" such configurations, but the 
number of independent probabilities is < 2", due to normalization, and various 
symmetries.) The system is assumed homogeneous, so that the P^"'^ are independent 
of position. 

Since transitions in a set of n contiguous sites generally depend on sites out- 
side the cluster, the n-site probabilities are coupled to those for n + 1 and so on, 
generating an infinite hierarchy. The n-site approximation truncates this hierarchy 
by approximating m-site probabilities (for m > n) in terms of n-site conditional 
probabilitites. For example, if fjfc represents the state of site k, then in the two-site 
approximation a three-site joint probability for a sequence of nearest-neighbor sites, 
k, I, m, is approximated so: 

P{ak, ai, am) ^ P{ak\ai)P{ai, am) = . (1) 

P{'^i) 

It is convenient to denote configurations in the PCP by a string of O's and I's, 
the former representing vacant sites and the latter, occupied sites. We denote the 
probability of configuration ai, ...,an by {ai, an), i.e., (1) denotes the probability 
of a randomly chosen site being occupied, (11) the probability of a nearest-neighbor 
occupied pair, and so on. 

In the one-site approximation we have the transitions (1) — >■ (0) at rate 2p(l)^ 
and (0) (1) at rate (1— p)(l)^[l — (l)]-|-/i[l — (1)]^, yielding the equation of motion: 
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-(l) = -2p(l)2 + (l-p)(0)(l)^ + M0f 



(2) 



where of course (0) = 1 — (1). For h = the stationary active solution is (1) = 
(1 — 3p)/{l — p) so that the critical annihilation probability Pc = 1/3. [Aside from 
the active solution, we always have the absorbing state, (1) = 0.] For h > 0, there is 
a nonzero stationary particle density for any p e [0, 1], i.e., the source removes the 
phase transition, at this level of approximation. 

Next we consider the pair approximation. There are two independent proba- 
bilities, since, by symmetry, (10) = (01), while normalization implies (00) + 2(01) 
+ (11) = 1. There are five possible transitions of a NN pair of sites: from (00) 
to (01); from (01) to (00) or (11); from (11) to (00) or (01). To illustrate how 
rates are calculated, we consider the transition (00) — > (01). There are two possible 
mechanisms: one involves a (11) pair just to the right of the central pair; the other 
involves the action of the source, and requires a vacant site to the right of the cen- 
tral pair. In the first case, the probability of the required configuration, (0011), is 
approximated as (00)(01)(11)/[(0)(1)]; the intrinsic rate is (1— p)/2. In the second 
case the configuration probability is (000) ~ (00)^/(0) in the pair approximation. 
Thus the contribution to d{01)/dt due to the transition (00) — > (01) is 

l-p(00)(01)(ll) (00)^ 
2 (0)(1) ^ (0) 

Note that this expression is multiplied by -2 in the equation for d{00)/dt, to take 
the mirror- image transition [(00) — > (10)] into account. (Recall that (10) has been 
eliminated by symmetry.) 

Proceeding in this manner, one readily obtains a pair of equations for (11) and 
(1) = (11) + (01): 

|(ii) ^ -mn) + mm + (i-P)ii - (11)] '""^-/"^' (3) 



l(i) = -.Kn).(i-,MM..M! 

The stationary active solution for h = is (1) = (1— 5p)/(l— p), (11) = (1— 3p)(l)/(l— p), 
so that Pc = 1/5. The critical annihilation rate is 2/3 for any nonzero h. 

The three-site approximation involves five independent variables and thirteen 
distinct transitions. Integrating the coupled equations numerically, one finds 
Pc — 0.1277 for h — 0, and Pc — 0.18197 for any nonzero source. We have ana- 
lyzed cluster approximations for up to 6 sites; the predictions for pc and for the 
density of isolated particles (pnat at the critical point are summarized in Table I. The 
cluster approximations appear to approach the simulation values (note the oscilla- 
tory nature of the approach to 0not)- We find in the n = 2 approximation that for 
nonzero h, the order parameter p oc {pc — pY in the neighborhood of the critical 
point, that is, the mean-field exponent (5mf = 2 in this case; for n > 2 we find the 
usual value, (3mf — 1- 
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For h — we can obtain a rongh estimate of Pc by extrapolating the cluster 
results, i.e., via a linear fit to Pc{n) plotted versus 1/n. Using the data for n = 
3 — 6 yields Pc = 0.082(7); using only the n — 4 — 6 data, we find Pc = 0.07(1), 
consistent with the simulation result. Unfortunately the behavior of these low-order 
approximations is not sufficiently regular to allow an extrapolation to n — > oo, when 
h>0. 

Given a sequence of cluster approximations, and knowing the value of pc, one 
may apply Suzuki's coherent anomaly analysis to extract certain critical exponents 
[43] . In the present case, however, our estimates for (3 appear to converge slowly; we 
find P ~ 0.29 ior h — (from a quadratic fit to the order parameter at the critical 
point), and ~ 0.28 for h = 0.5. Data for larger clusters will be needed in order to 
derive precise predictions for critical exponents. 

In Fig. 1 we compare the n-site approximation predictions against simulation 
results for the phase boundary, Pc{h). The cluster approximations appear to ap- 
proach the simulation curve in a qualitative fashion. On the other hand, all the 
approximations (n < 6) predict a jump in pc at h = 0, while simulations show pc to 
be continuous, though singular, at this point. Thus it appears that detailed features 
of the phase diagram are beyond the small-cell approximations developed here. 

IV. SIMULATIONS 
A. Method 

We first define the simulation algorithm for the PGP in the absence of a source 
{h = 0). We consider a ring of L sites. Since all events depend upon the presence 
of a NN pair, we maintain a list of such pairs. An event consists in (1) choosing a 
process (annihilation with probability p, creation with probability 1— p); (2) choosing 
a NN pair at random from the list; (3) deleting the two particles in the pair, in case 
of annihilation, or, in case of creation, adding a new particle adjacent to the pair, 
if possible. More specifically, in the case of creation, step (3) consists in choosing 
(with equal probabilitites) the site immediately to the left, or to the right, of the 
NN pair, and inserting a particle at this site if it is vacant. The time increment 
associated with each event (creation or annihilation) is At — 1/Np where Np is the 
number of NN pairs just prior to the event. (If A'p = the system has fallen into an 
absorbing configuration and the trial ends.) In this way each NN pair undergoes, 
on average, one event (creation or annihilation) per unit time. The list of NN pairs 
must naturally be updated following each annihilation or sucessful creation. 

When h > 0, each creation or annihilation is followed by a certain number Nh 
of source events. Suppose a creation or annihilation has just occured, and that 
there are now Np pairs. Then the next creation or annihilation will occur after a 
time interval of l/Np] during this interval, the expected number of source events is 
LhAt — Lh/Np, since h is the rate of insertion attempts per site and unit time. Thus 
after each creation or annihilation we perform Nh insertion attempts, where Nh is 
a Poissonian random variable with mean Lh/Np. An insertion attempt consists in 
choosing a site i at random (each site has probability 1/L; no lists are used in this 
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process), and inserting a particle at there if and only if sites i — 1, i, and i + 1 are 
all vacant. 

B. Results 

We used the simulation algorithm described above to determine Pc{h) at a num- 
ber of different h values. In addition, detailed studies were performed ai, h — 0, 0.1, 
and 10 in order to obtain precise values for critical exponents and other properties. 
All simulations start from a fully occupied lattice. The order parameter p = Np/L 
(i.e., the density of nearest-neighbor pairs), the survival probability, P{t) (that is, 
to have Np > 0), and the particle density, ip, are monitored during the evolution, 
and, in particular, in the stationary state. The properties of principal interest as 
regards critical behavior are the stationary order parameter p (the overlinc denotes 
a stationary mean), the moment ratio m = p^/ij^Y, and the hfetime r, obtained 
by fitting an exponential to the survival probability: P(t) ~ e^*/^. Uncertainty 
estimates for p, m and r are obtained from the standard deviation evaluated over 
a set of 3-5 independent runs. 

In order to fix the critical point, foi' a given value of h, we employ three 
criteria. The first is the finite-size scaling behavior of the order parameter: one 
expects a simple power-law dependence at the critical point, 

p(Pe, L) ~ L-P'^- , (5) 

as has in fact been verified for many absorbing-state phase transitions [2,3]. For 
noncritical p-values, the dependence of p(pc; L,) upon L shows deviations from a 
power law, typically manifested as curvature on a plot of p(pc) L) versus L on log 
scales. The same considerations apply to the survival time, whose critical finite-size 
scaling form is 

T{p,,L)^U\\'^^ . (6) 

Our second criterion for locating p^ is thus power-law dependence of r on the system 
size. 

A third criterion is based on the independence of moment ratios such as m at 
the critical point. This property, which has long been used in studies of equilibrium 
critical phenomena [44], was more recently verified for absorbing-state phase tran- 
sitions in the contact process and PGP [32,45]. We determine the value p*{Li,L2) 
at which the moment ratios m(p, Li) and m{p, L2) take the same value, for pairs of 
system sizes Li and L2. Extrapolating p* to L — >■ 00 yields an estimate for p^.. 

We studied system sizes L — 50, 100, 200, 500, and 1000. A study at a particular 
value of h begins with a quick survey of small systems, to obtain a preliminary 
estimate of Pc- Then, for each system size, we perform high-statistics studies at 
three to five p- values, obtaining the properties of interest to high precision. Values 
of p, m, r, and ^ at intermediate p-values are obtained via interpolation of the 
high-precision data using polynomial (typically cubic) least-squares fits. 
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We tested this procedure on the standard PCP {h = 0). For L = 50 we ran a total 
of 6 X 10^ trials at p-values of 0.0769, 0.0770, 0.0771, and 0.0772. For this small 
system the stationary state is well sampled in trials that extend to a maximum time 
of 500. For L= 1000, we performed a total of 3 x 10^ trials (extending to a maximum 
time of 10^) at p = 0.0770, 0.0771, and 0.0772. Sample sizes and maximum times 
for intermediate system sizes fall between the values quoted for L = 50 and 1000; 
similar parameters were used in the high-precision studies at h = 0.1 and 10. Figure 
2 illustrates the results of the fitting procedure for p; it is clear that the simulation 
data are well represented by a smooth function. In Fig. 3 we show all of the data 
for m (for h — 0) together with the associated polynomial fits. 

For h = 0, the three criteria mentioned above yield pc = 0.077091(5), in excellent 
agreement with previous studies [32,40,46]. The uncertainty in pc reflects scatter 
in the last three moment-crossing points, and in the p-value yielding power laws 
for p and for r. Extrapolating the values of m{pc, L) to L ^ oo, we find m = 
1.1740(5). (Here the uncertainty includes three sources: that in the original data, 
the uncertainty of the extrapolation at pc, and the uncertainty in pc itself. The 
same applies to the exponent values discussed below.) Our result is once again in 
excellent agreement with that of Ref. [32], rric = 1.1735(5) for transitions in the DP 
universality class in 1-1-1 dimension. 

From the scahng of p at the critical point we find P/iy± — 0.2522(5) (see Fig. 4); 
this exponent ratio takes the vahic 0.2521 for the DP class [47]. The survival time 
data (Fig. 5) yield = 1.577(4), while the DP value is 1.5807 [47]. As noted 

in the Introduction, the order parameter is expected to follow a power law, p ~ t~^, 
during the approach to the stationary state. A study at p = 0.077091 in a system 
of 1000 sites yielded 9 — 0.1596(2) (see Fig. 6), while the expected DP value is 
9 = j3/v\\ = 0.15947(3). Finally, a study of p in the supercritical regime (p = 0.0722 
- 0.076, for L = 1000) yielded p = 0.276(3) [via the usual relation, p ^ {pc- pf]- 
The accepted value is /? = 0.27649(4) [47]. We note that these exponent values 
are the highest precision yet reported for the PCP, and that they place the model 
unambiguously in the DP universality class. 

While the results for /S/z/j^ and P furnish the estimate z/^ = 1.092(7), it is de- 
sirable, for reasons explained below, to have an independent estimate of z/_l. Since 
A = Pc — p enters all finite-size scaling forms in the combination AL^^"^ , wc can 
obtain such an estimate either from a data-collapse analysis, or by studying the 
finite-size dependence of a derivative such as dm/ dp at Pc] the latter must diverge as 
. "\Yc find that this derivative does follow a power law; a fit to the data yields 
z/_L = 1.086(22), where the relatively large error refiects the uncertainty associated 
with numerical evaluation of a derivative. (Both our estimates are consistent with 
the DP value, = 1.0968.) _ 

Our data for the particle density V(Pc, L) fall on a straight line when plotted 
versus L"^/"-^, and extrapolate to 0.241(1) for L— >oo (Fig. 7). In the active state, 
i/j represents the total particle density (isolated particles as well as those belonging 
to pairs), but when L^oo at p = Pc there are no pairs; in this limit ip corresponds 
to the natural density, 4>nat, of isolated particles. Our result agrees with previous 
studies, based directly on absorbing configurations [40,41]. 
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The procedure described above, which furnishes results of unprecedented preci- 
sion for the standard PCP, was repeated for h = 0.1 and 10. The resuhing critical 
parameters are compared with the h = case in Table II. We remark that the depen- 
dence of p, T, and m on L, at the critical point, is qualitatively the same as for h — 0, 
and that the associated critical exponents, /3/i'±, i'\\/i'±, and 9, and moment ratio 
m, are the same, to within uncertainty, as for h = 0. This is particularly clear from 
Fig. 4, which compares the dependence of p on L at the critical point for h = 0, 0.1, 
and 10, and Fig. 5, a similar plot of the lifetime r. The same power laws are seen, 
regardless of the value of h, the sole difference between zero and nonzero h being a 
slight change in amplitude. The data for h — 0.1 and h — 10 are virtually identical. 

The decay of p at the critical point again appears to follow the same power law, 
regardless of whether or not the source is present. Fig. 6 shows that the main effect 
of the source is a small decrease in the amplitude. Analysis of the local slopes of 
the graphs in Fig. 6 leads to ^ = 0.1603(5) in the presence of the source, in very 
good agreement with the value 9 — 0.1596(4) found for h — 0. The rather imprecise 
estimates for u^, obtained directly from the scaling of dm/ dp, likewise agree to 
within uncertainty, independent of h. 

There are two principal differences between the scaling properties with and with- 
out the source. First, Table II shows that the hmiting (L— >oo) value of the particle 
density ip is much larger when h > 0. (The data for ip at the critical point, and the 
extrapolation to ipoo, are depicted in Fig. 7.) In fact, there is good reason to expect 
V'oo to jump from (pnat to a higher value as soon as h > (just as in the cluster 
approximations of the previous section), since, no matter how weak, the source will 
eventually fill in all available sites (as in RSA), in the absence of activity. With 
increasing h, V'oo approaches the RSA value of 0.4323. 

In addition to the expected difference in limiting particle densitites, the source 
appears to induce a rather subtle change in the critical exponent (3 (see Fig. 8). 
The value for /i > is about 4% greater than for h = 0; since only about 2% of 
this difference can be attributed to uncertainty, this seems to represent a significant, 
albeit small, change in the exponent. The interpretion this result, however, is not 
straightforward, since all of the other exponents (or ratios) studied appear to be 
insensitive to the source. In particular, the values for jS/u^ are constant to within 
1.2%, while our result for i>± (from the scaling of dm/ dp) is constant to within 0.6%, 
implying constancy of /3 to better than 2%. Thus it is difficult to accept the apparent 
change in /3 at face value; a better understanding will require either a theory of the 
effect of the source, or numerical results for larger systems. It is worth noting that 
with the exception of f3, all exponents have been obtained via finite-size scaling at 
the critical point. Thus our results are consistent with the possibility that the source 
modifies the apparent scaling of the order parameter near, but not at, the critical 
point, perhaps due to additional corrections to scaling. We cannot rule out a true (as 
opposed to merely apparent) change in the exponent /3 with the data in hand, but, 
given our other results, such a modification would entail a very surprising violation 
of the usual connection between scaling in the supercritical regime and finite-size 
scaling at the critical point. 

While the effect of the source on critical exponents, if any, is quite subtle, its 
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effect on tfie position of tlie critical point, Pc, is dramatic. Fig. 9 shows the phase 
boundary in the h-p plane, as determined in simulations. Evidently Pc grows in 
a singular manner as h is increased from zero. Although there is no jump in p^, 
(contrary to the predictions of the cluster approximations described in the preceding 
section), all evidence points to the derivative dp/dh being infinite at /i = 0. For 
larger h, Pc saturates; plotting pc versus 1/h leads to the estimate Pc = 0.0966(1) in 
the h^oo limit. The nature of the singularity at h = is not fully clear: a double 
logarithmic plot of Apc = Pc{h) — Pc(0) versus h does not yield a simple power law, 
perhaps due to saturation-induced curvature, but suggests an exponent of ~ 1/2 
as /i — > (see Fig. 9, inset). As can be seen from Fig. 1, it is plausible that Ap 
is a smooth function of h^/"^. Elucidating the precise nature of the singularity is a 
central issue for future study. 

V. DISCUSSION 

Motivated by the question of universality at an absorbing-state phase transition 
in the presence of a source coupling to a non-ordering field, we investigated the 
phase diagram and critical behavior of the PCP with a source of isolated particles. 
We studied the model via cluster approximations and extensive simulations, n-site 
cluster approximations (for n < 6) yield predictions for the phase boundary Pc{h), 
and the isolated particle density (p, that appear to approach the simulation results. 
All of the cluster approximations studied here predict of a discontinuity in pc at 
h = 0; simulations show Pc{h) to be a continuous function, albeit singlar ai h — 0. 

The central conclusion from extensive simulations performed for source strengths 
/i = 0, 0.1, and 10 is that the presence of the source, and the associated change in 
the background density of isolated particles, has no detectable effect on scaling at 
the critical point. The present level of precision permits us to state that the values 
of i'±, l3/i'±, Q-i and the moment ratio m — fp'/ifif^ are constant to within 

1% as we vary the intensity of the source. Our results for these parameters are fully 
consistent with known values for the DP universality class. We have noted, on the 
other hand, that in the presence of the source, the critical exponent /3 appears to 
be about 4% greater than the DP value. This is in conflict with the balance of our 
results which indicate constancy of the critical exponents. A resolution must await 
development of a theoretical understanding of scaling in the presence of the source, 
and/or more extensive simulations. 

The apparent insensitivity of the critical behavior to a change in the stationary 
density places the PCP in the same category as the threshold transfer process (TTP) 
[27,38]. Whereas in the TTP the nonordering field can relax in the absence of activ- 
ity, while in the PCP it cannot, the two models are similar in that the nonordering 
field exhibits only short-range correlations [26,29]. An interesting open question is 
whether perturbations in a nonordering field that lead to long-range correlations can 
affect critical behavior. 

Our simulations reveal that the phase boundary, Pc{h), is singular at h — 0. The 
available data are consistent with an asymptotic power law Ap^ ~ h" with a ~ 1/2, 
but further studies will be needed to characterize the singularity precisely. A scaling 
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or renormalization-group analysis of the PCP with an isolated particle source would 
clearly be desirable, in order to understand the form of the phase diagram, and the 
scaling of the order parameter. 
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TABLES 



TABLE I. n-site approximation results 



n 


Pc{h = 0) 


Pc{h = 0.1) 


^{h = 0) 


(/)(/i = 0.1) 


2 


0.2000 


0.6667 





0.5 


3 


0.1277 


0.1820 


0.23 


0.461 


4 


0.1185 


0.15023 


0.2166 


0.5 


5 


0.114 


0.14405 


0.18 


0.46 


6 


0.103 


0.117 


0.2234 


0.438 


SIM 


0.077091(5) 


0.086272(15) 


0.241(1) 


0.421(1) 



TABLE XL Critical parameters of DP and the PCP. DP exponents from Ref. [47], m from 
Ref. [32]. Numbers in parentheses denote uncertainties in the last figure(s). 



h 


Pc 


f3/i^± 


1^1 |M 


f3 




e 


m 


-000 


DP 




0.25208(4) 


1.5807 


0.27649 


1.09684(1) 


0.15947(3) 


1.1735(5) 






0.1 
10 


0.077091(5) 
0.086272(15) 
0.097850(10) 


0.2523(3) 
0.2540(15) 
0.2554(10) 


1.577(4) 
1.571(9) 
1.574(10) 


0.276(3) 
0.287(3) 
0.287(3) 


1.086(22) 
1.083(20) 
1.080(10) 


0.1596(4) 
0.1602(4) 
0.1604(4) 


1.1740(5) 
1.176(1) 
1.175(1) 


0.241(1) 
0.421(1) 
0.433(1) 
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Figure Captions 



Fig. 1. Cluster approximation predictions (for n = 4, 5, and 6, as indicated) and 
simulation results (circles) for the phase boundary. 

Fig. 2. Simulation data for the stationary order parameter value, p, for L = 50 and 
h = 0. The solid line is a cubic fit. 

Fig. 3. Simulation data and polynomial fits for the moment ratio m, for h = 0. 
System sizes L = 50, 100, 200, 500 and 1000 in order of increasing slope. Inset: 
detail of crossing region for L = 100,. ..1000. 

Fig. 4. Stationary density versus system size at the critical point for h — (□), 
h — 0.1 (+), and h — 10 {<>). Error bars are smaller than the symbols. 

Fig. 5. Lifetime versus system size at the critical point. Symbols as in Fig. 4. 

Fig. 6. Decay of the order parameter at Pc- Upper curve: h = 0; dotted line h = 0.1; 
lowermost curves: h = 1 and h = 10, which cannot be distinguished on this scale 

Fig. 7. Stationary particle density versus L'^^'^-^, illustrating linear extrapolation 
to -000- Symbols as in Fig. 4. 

Fig. 8. Stationary pair density versus A — Pc — p, ior h — 0, 0.1, and 10; L = 1000. 
Symbols as in Fig. 4. Straight lines represnet least-squares fits to the linear portion 
of the data ior h — and h — 10; the slopes are 0.276 and 0.287, respectively. 

Fig. 9. Critical annihlation rate Pc versus source strength h (the dashed line is 
simply a guide to the eye). The inset shows the critical point shift Apc versus h on 
log scales; the straight line has a slope of 1/2. 
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